Decay Rate Distributions of Disordered Slabs and Application to Random Lasers 
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We compute the distribution of the decay rates (also referred to as residues) of the eigenstates of a 
disordered slab from a numerical model. From the results of the numerical simulations, we are able 
to find simple analytical formulae that describe those results well. This is possible for samples both 
in the diffusive and in the localised regime. As example of a possible application, we investigate the 
lasing threshold of random lasers. 



I. INTRODUCTION 

A very successful approach to describe disordered ma- 
terials is supplied by random-matrix theory, see Refs. [|,|2| 
for reviews. While one can put the beginning of random- 
matrix theory at Wigner's surmise for describing the 
scattering spectra of heavy atomic nucleijl its theoret- 
ical foundations were laid only much later .□ It was very 
successfully applied to electronic transport in disordered 
wires and mesoscopic quantum dots, and recently these 
methods have been adopted to model (quantum) trans- 
port of optical radiationJa|jaedia with spatially fluctu- 
ating dielectric constant SBOu 

In the theoretical treatment of disordered materials, 
two particular geometries are of special importance, 
namely the disordered slab and the chaotic cavity (see 
Fig. |l|). The principal difference between the geome- 
tries is easily explained: A chaotic cavity is an object 
in which the dynamics is chaotic due to the shape of the 
cavity or due to scatterers placed at random positions. 
The size of the opening is small compared to the total 
surface area of the cavity. Particles (electrons, photons) 
are then "trapped" inside the structure for a time that 
is long enough to crgodically explore the entire cavity. 
In a disordered slab, particles cannot be trapped that ef- 
ficiently. They can no longer explore the entire volume 
ergodically but they still stay long enough to explore the 
entire cross-section of the sample, thus still making a 
random-matrix description possible. In order to call this 
geometry a "wire" or a "waveguide" its length should 
be much larger than its width. To be able to apply the 
theory only the much weaker criterion that the length is 





FIG. 1: Two frequent geometries in the theory of disordered 
media are the chaotic cavity (left) and the disordered slab 
(right). The motion in the chaotic cavity is completely ergodic 
while in the disordered slab it is ergodic only in the transversal 
direction. 



sufficiently larger than the mean-free path of the medium 
has to be fulfilled. 

Two different aspects are of special important in the 
theory of disordered media, namely transport properties 
and resonances. The transport properties (moments of 
the eigenvalues of the transmission and reflection matri- 
ces) are-known for the disordered slab in the limit that it 
is wideja for the chaotic cavity with an opening that is so 
snaali-that only one or two modes can propagate-through 
itJiO'Eil or a chaotic cavity with a wide openingE 

Less is known about the poles of such systems. (The 
eigenvalues of the Hamiltonian correspond to poles of the 
scattering matrix, and these show up as resonances in a 
scattering experiment. Hence the somewhat inconsistent 
nomenclature.) The beginning of random-matrix theory 
can be put at the moment when Wigner surmised,rth6| 
eigenvalue distribution for a closed chaotic cavitylffl'tHl 
Here we are interested in open systems, where the eigen- 
values acquire an imaginary part. (The imaginary part 
is referred to as residue.) It determines the decay rate 
of the (quasi-)eigenstate of the system. For chaotic cavi- 
ties with broken time-reversal symmetry, the decay rate 
distribution, is known analytically for an opening of arbi- 
trary size.Ej The distribution for the more important case 
of preserved time-reversal symmetryL^I is not known but 
can be approximated by a cavity with broken symmetry 
and an opening of half the real size. 

Information on the residues of a disordered slab is 
very limited, and only the scaling behaviour of the 
large residue-tail in the localised regime was determined 
recentlyJiju3 This deficiency is felt especially strong in 
the random-laser community since the location of the 
residues gives the lasing threshold of an optical system, 
and most experimental setups resemble a disordered slab 
much more than they resemble a cavity. This paper fills 
this gap by presenting the results of numerical simula- 
tions. The quality of the numerical decay rate distribu- 
tions is good enough that it allows us to arrive at analyt- 
ical formulae for the distribution function, including its 
dependence on the parameters of the system. The idea 
to use high-quality simulations to arrive at formulae is 
not completely new as the distribution of the scattering 
strengths of chaotic cavities was found in the same way.El 

This paper is organised as follows: First, we introduce 
the Anderson Hamiltonian used the describe the disor- 
dered slab. In Sec. Ill we show how the eigenvalues of 
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that Hamiltonian can be computed in a efficient numer- 
ical way. Depending on the length of the slab, it can be 
in either the diffusive or in the localised regime. We will 
analyse the decay rate distributions for both regimes sep- 
arately, first in Sec. ^ for the diffusive and then in Sec. V. 
for the localised regime. Until that point all results are 
completely general and ca n be applied to electronic and 
photonic systems. In Sec. VII we specialise on the las- 
ing threshold in amplifying disordered media. We dis- 
tinguish between the diffusive and the localised regimes 
(Sees. [VlTA] and[viTB|). We conclude in Sec. |V1II. 




II. ANDERSON HAMILTONIAN FOR A 
DISORDERED SLAB 

We consider a two-dimensional slab of length L and 
width N. The slab is described by an Anderson-type lat- 
tice Hamiltonian with spacing I. In the Anderson model, 
transport is modelled by nearest-neighbour hopping be- 
tween lattice sites. Without loss of generality we can set 
the hopping rate to I. With a spatially varying potential 
P(x,y) the Hamiltonian becomes 





= P(x,y) 
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All other elements are zero, x runs from 1 to N, and y 
from 1 to L. 

The real part E of the eigenvalues of Ti in the limit 
of large N and L is confined to the interval [—4; 4]. (If 
the average of P(x, y) is nonzero, the interval is simply 
shifted by that average. If P(x, y) is fluctuating as a func- 
tion of x and y — like it does for a disordered medium 
— the interval becomes a bit wider.) For electronic sys- 
tems, E gives the energy of the eigenstate, and Eq. (|l|) 
hence describes a slab with a conduction band of width 
8. For photonic systems, the real part of the eigenvalue 
gives the eigenfrequency. For both systems, the imagi- 
nary part 7 of the eigenvalue gives the decay rate of the 
eigenstate. (Actually not 7 but rather 7/2 is the decay 
rate but for the ease of notation we will continue to refer 
to 7 simply as the decay rate.) 

k in Eq. (|lb|) quantifies the strjength of the coupling 
between the slab and the outsideJl3 Using the units in- 
troduced above, k has the value sin 2 k where k is the 
wavevector at the energy at which particles are injected 
respectively emitted. This quantity is proportional to the 
velocity of the particle perpendicular to the interface. In 
the centre of the band sinfc = 1 whereas at the edges 
sin/c = 0. 

If k is chosen to be constant (i. e. not to depend on 
energy) ideal outcoupling can be described only for one 
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FIG. 2: The Hamiltonian TC has a band structure. The thin 
lines contain matrix elements that are mostly 1, the diago- 
nal is filled with complex numbers, and all other elements 
are zero. The entire matrix is symmetric but non-Hermitian 
(since there are complex elements on the diagonal). 



specific value of the energy. We will do this since oth- 
erwise solving the Hamiltonian no longer is a standard 
eigenvalue problem, and set k = I .-hence modelling ideal 
coupling at the centre of the band£3 Working at the cen- 
tre of the conduction bands offers the advantage that 
the width N of the sample is identical to the number N 
of propagating modes, and thus allows the describe the 
largest number of propagating modes for given size of 
the Hamiltonian (i. e. given numerical work). Itis possi- 
ble to include energy dependent coupling termsEj but it 
should be stressed that a constant k is more efficient and 
gives completely correct results as long as only eigenval- 
ues near the respective energy are considered. We set 
k = 1, meaning ideal coupling at the centre of the con- 
duction band. 

It should be stressed that - even though we are mod- 
elling a two-dimensional system — the results are valid 
for three-dimensional systems as long as L > N. A par- 
ticle that is injected into the slab ergodically explores 
the entire cross-section of the sample before being emit- 
ted again, and hence loose its memory of which sites are 
"connected" by hopping terms. The sites can then be 
re-arranged, e. g. in a three-dimensional structure. Only 
for very short samples, L < N , this is not possible but 
for such samples already applying the Anderson model 
(i.e. only allowing nearest-neighbour hopping) is very 
questionable. The only "real" restriction that can limit 
the application of our results to certain photonic three- 
dimensional systems is that particles can leave the sample 
only at the front and at the back — and not through the 
"sides" . 

In the formulation of Eq. ([j]) the matrix H. has double 
indices but these are easily removed by considering Ti nn ' 
with n = x + (y — 1)A^. (It would not make sense to 
set n — y + (x — 1)L since usually L » N, and we 
want to work with a band matrix that is as small as 
possible.) This results is a matrix of the form as depicted 
in Fig. 0. It is a banded LN x LN matrix with band 
width 2N + 1. Also within the band most elements are 
zero (since usually N 3> 1). The matrix is symmetric 
but non-Hermitian as there are complex numbers on the 
diagonal. 

The model ([!]) has been widely used since an efficient 
way to compute the transmission through such a slab 
is known. E3 The method of recursive Greens functions 
allows to compute the entire scattering matrix, hence all 
linear transport properties, in a time of order 0(LN 3 ) 
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and with only minimal storage requirements 0(N 2 ). No 
explicit reference to the Hamiltonian Ti. is made, so that 
eigenvalues cannot be computed by this method. 

III. COMPUTING EIGENVALUES OF 
SYMMETRIC COMPLEX NON-HERMITIAN 
BANDED MATRICES 

Since the Hamiltonian H from Eq. ([j]) is both banded 
and sparse one might be tempted to use an eigensolver 
for sparse matrices to compute the eigenvalues of Eq. (Q) . 
A sparse eigenvalue routine needs to be able to solve the 
equation 

(H-fil)x>=y (2) 

for the unknown vector x for arbitrary /i and y. In partic- 
ular, the eigensolver needs to set \x close to an eigenvalue 
of H so that the matrix Ti — /xl is ill-conditioned. A 
numerical solution of Eq. (|J) is then difficult and very 
expensive. Furthermore, only one eigenvalue is found 
at a time, and control of which eigenvalue the algorithm 
will converge to is difficult. (Algorithms for sparse matri- 
ces always use inverse iteration so that the corresponding 
eigenvector will be returned without additional effort but 
the eigenvector is of no use for us.) Using an algorithm 
for banded matrices is thus the better alternative. 

There exist a number of algorithms for real symmetric 
or complex Hcrmitian band matrices. Both problems are 
characterised by real eigenvalues, so that they are con- 
ceptually identical. Only one algorithm for computing an 
eigenvalue (plus the corresponding eigenvector ]_of a gen- 
eral complex band matrix has been published. E3 It uses 
inverse iteration, so it is of limited use here. 

We thus had to implement our own eigenvalue solver. 
The eigenrepresentation of TL in terms of the diagonal 
matrix A = diag(Ai, . . . , Xjq) of the eigenvalues \ of A 
and the matrix U of eigenvectors is 

U = uau- 1 . (3) 

We now observe that for symmetric, that includes com- 
plex symmetric, 7i it is always possible to chose eigen- 
vectors such that UU T = 1. If U would be a real matrix, 
one would call U orthogonal but since it is complex there 
is no special name for the property UU T = 1. 

Algorithms for diagonalising a real symmetric matrix 
A implicitly decompose A as 

A = QDQ T , QQ T = 1 (4) 

with the matrix D of eigenvalues of A. It is therefore 
possible to adapt such an algorithm for our needs. Most 
algorithms for banded matrices first reduce A to tridiago- 
nal form A' by transformations of the form A' — Q'AQ' T , 
and we will adopt this strategy. (A matrix is called tridi- 
agonal if only the diagonal and its neighbouring elements 
are nonzero. If A would be real, the transformation 
A — » A' would be called a similarity transformation.) For 



a band matrix this is possible in an efficient way since it 
is not necessary to compute (and thus store) the full ma- 
trices Q', and by annihilating the elements of the matrix 
A in a clever order, the band structure is kept intact in 
all steps Jl3'E3 

The reduction of the complex matrix H to tridiagonal 
shape is done by straight-forward adaptation of this al- 
gorithm from real to complex numbers, where care needs 
to be taken that the dot product x • y = J2i x iVi IS use d 
and not the dot product x ■ y = "Y^^Xiyi normally used 
for complex vectors. (The overbar marks the complex 
conjugate.) 

To compute the eigenvalues of the tridiagonal matrix 
for the real symmetric or complex Hermitian case, meth- 
ods that isolate eigenvalues in disjunct intervals, are used 
("divide and conquer" and similar methodsEil). Such 
methods work for both of these cases as all eigenvalues 
are real and can thus be ordered. This no longer is pos- 
sible here as the eigenvalues are conopJex. We therefore 
use the QR respectively QL method£3 

For an K x K banded matrix with band width W 
the time needed to compute the eigenvalues increases as 
0(K 2 W) whereas the storage requirements increase as 
0(KW). In terms of the dimensions L and N of the 
disordered slab, this means that the time increases as 
0(L 2 N 3 ) and the storage space as 0(LN 2 ). For given 
computational resources, both scalings impose an upper 
limit on the system size that can feasible be treated. For 
typical values of the ratio L and N, and for "realistic" 
computer equipment, the time .limit is reached somewhat 
earlier than the memory limit. l3 

With respect to a similar algorithm for full matrices 
one wins a factor L (usually of order several hundred) for 
both time and memory by using the banded algorithm, 
thus allowing to treat system that could not be treated 
otherwise. Still, the work presented in this paper is a big 
numerical challenge. To arrive at the results, of the order 
of 100000 hours of cpu time on fast PC's were needed. 



IV. NUMERICAL SIMULATIONS 

Disorder is modelled by assigning random values to 
P(x,y). It is assumed that those random numbers are 
uniformly distributed in the interval [—w; w] so that w 
measures the amount of disorder. 

We only consider eigenvalues near the centre of the 
conduction band as the assumption of ideal coupling is 
only valid there. For numerical reasons it is essential that 
the centre of the conduction band is at £. = 0, i.e. one is 
not allowed to add an offset to P(x, y)£^ 

We hence chose a window [—d; d] and only include 
eigenvalues in the further analyses when their real part 
is inside that window. If the window is too large, sys- 
tematic errors are introduced while too small a window 
leads to bad statistics. As can be seen from Fig. for of 
d = 0.1 the distribution function already agrees with the 
distribution function for d = 0.01 but has much better 
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statistics, d — 0.5 and d — 1.0 gives significant system- 
atic deviations. For this reason, all results presented in 
this paper assume a window with d = 0.1. 

The formulation of the model in Sec. || is in terms of 
generic units. Contact with a microscopic model or an 
experiment is best made in terms of the mean free path. 
It can be computed from the length-dependence of the 
transmission probability T through the sample. In the 
diffusive regime, I < L <C Nl, it is given b 



3" 



l 

f 



L 
1 



(5) 



The mean free path can be computed by fitting T(L) to 
this functional form. 

The transmission probability has been computed using 
the method of recursive Green's functional^ for variable 
disorder strength w. As Fig. || shows, the numerically 
computed mean free path I is for the range of w in ques- 
tion in very good approximation given by 



I = 



,,3/2 



(0) 



(Computed for each value of w from 50 samples with 
L = 2, 4, . . . , 98, 100 and N = 50.) In the following, we 
will no longer make explicit reference to w but rather 
give the more intuitive mean- free path I. 



V. DIFFUSIVE REGIME 

For a sample length L with I < L <C Nl the sample 
is said to be in the diffusive regime. It is immediately 
obvious that the diffusive regime can only be observed in 
sufficiently wide samples, N 3> 1. 

For chaotic cavities with broken time-reversal symme- 
try an analytical result for the decay rate distribution has 
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FIG. 3: Probability distribution of the decay rates for given 
system parameters as a function of the window size around 
the centre in which eigenvalues are included in computing the 
probability distribution. The different dots mark the distri- 
butions with d = 0.01, 0.1, 0.5 and 1.0. 
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FIG. 4: Mean free path Z as a function of disorder strength w 
from a numerical simulation (crosses) and from Eq. (^). 

been given by Fyodorov and Sommersill We start from 
their result and rescale it, 



y 2 Ml 



My 



dx = 



y 



M 



(7) 

V[y) is normalised to one and in our scaling is for all 
M > 1 peaked near a value of y of order 1. In the origi- 
nal formulation for a chaotic cavity, M is the number of 
modes propagating through the opening of the cavity. 

In the following we will argue that the decay rate dis- 
tribution P{g) can be written in the form (m) as 



Pin) = -v(^: 

7o 7o 



(8) 



with some scaling factor 70 and some effective number 
of modes M ^ N. In Fig. | a comparison between the 
analytical suggestion and a simulation is given, and the 
agreement is striking. The horizontal axis has been plot- 
ted logarithmically since this results in both the differ- 
ences between the P(y) for different N becoming easier 
to recognise and in giving a more prominent place to the 
small-7 tail of P(7). In most applications, including the 
random laser discussed later in this paper, one is much 
more interested in small 7 than in large 7. 

The results of the simulations are fitted "by eye" 
against the functional form (Q), resulting in one pair of 
values for 70 and M for each set of parameters. Espe- 
cially at very small 7, there are sometimes numerical er- 
rors that introduce artifacts into the numerical histogram 
so that using an automatic fitting algorithm is not fea- 
sible. (Usually we computed 500-1500 realisations for 
each parameter set.) From our simulations, we find that 
the scaling factor 70 only depends on the length L of the 
sample and its mean free path I but not on its width N, 
and seems to be given by 



70 



21 
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(9) 
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FIG. 5: The numerically computed probability distribution 
P(g) for L — 175, N = 50, / = 12.9 and comparison with 
Eq. (Q) with M = 16. 



As Fig. g shows, the agreement between the result of the 
numerical simulations and Eq. (|) is good, and all ma- 
jor deviations are for small L where universal scaling is 
expected to be worse than for larger L. The model equa- 
tions set the speed of propagation to 1 but it is obvious 
that for some other choice for the propagation speed c 
one has to change Eq. (^) to 70 = 2cl/L 2 . 

While the determination of 70 is very precise, there 
is a somewhat larger error involved in determining M by 
fitting the analytical form to the results of numerical sim- 
ulations. First, we only fitted against integer M, though 
in principle a generalisation of Eq. (R) to noninteger M 
is possible, see Eq. (25). Secondly, if M > 25, the dif- 
ference between P(y) for M and for M + 1 becomes too 
small to tell with certainty which of these two values de- 
scribes the numerical result better. Thirdly and finally, 
even with 500-1500 samples for each set of parameter val- 
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FIG. 6: Scaling factor 70 as a function of I and L for Z = 24 
and N = 10, 20, 50, 100 (open symbols) and I = 12.9 and 
N = 50 (solid circles). The line marks the prediction from 
Eq. EL 



FIG. 7: Effective number M of degrees of freedom as a func- 
tion of the length L of the sample for a sample with N = 30, 
N = 50 respectively N — 70 propagating modes. The line 
marks the prediction from Eq. hd, the points are from numer- 
ical simulations. The size of the "errorbars" does not indicate 
some estimated error interval but simply marks the computed 
value ±1. 



ues, there are still some fluctuations in the numerically 
computed histogram for the decay rate distribution that 
in some cases make the decision on the right M a bit 
difficult. Considering all of this, one should allow for an 
error of 1 for M, and even of 2 for M > 25. 

We have computed M for a series of samples with in- 
creasing length for three different widths N. As Fig. |?j 
shows, the effective number M of modes is well approxi- 
mated by 



M = 



N 



l + L/{U) 



(10) 



The agreement between this suggested analytical form 
and the numerical simulation becomes better as the 
width N of the sample is increased. From the simula- 
tions it is obvious that the functional form Eq. (|h]) is 
correct but there still is the (small) possibility that the 
factor 6 might need to be replaced by a slightly smaller 
value. To answer this question with certainty, we would 
need to increase both L and iV significantly. Unfortu- 
nately, such simulations are outside the present time and 
memory constraints. 

Equations (Q-|l(J) give a good description of the de- 
cay rate distribution of a disordered slab in the diffusive 
regime, provided the slab is sufficiently wide. Since the 
transversal length scales are set by microscopic quantities 
(wave length of the light for optical systems, Fermi wave 
length for electronic systems), all macroscopic objects are 
"wide". 



VI. LOCALISED REGIME 

If the length L of a disorder medium is increased, the 
phenomenon of localisation sets in once L > Nl (see 
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FIG. 8: Numerically computed distribution of the decay rate 
7 for a sample in the localised regime (L = 71.55 I, N — 15) 
and comparison with the log-normal distribution Q13|) . 

Rcf. [j] for a review). In the localised regime the prob- 
ability of transmission T through the sample is reduced 
significantly and decays exponentially with the length L 
of the sample. The length scale £ is called the localisa- 
tion length, and can be computed from an ensemble of 
disordered slabs by computing the average of the loga- 
rithm of the transmission as a function of the length of 
the samples, hence 



L/£=QnT(L)) 



(11) 



One should note that this is not identical to fitting the 
transmission to (T(L)) oc exp(— L/£) since the large 
sample-to-sample fluctuations of T(L) in the localised 
regime would give a value for £ that is off by a factor 
4. The localisation length can also be computed analyti- 
cally from the Bafiafojy^ pfctfa using the DPMK equation, 
with the resulted 



N+l 



(12) 



and agrees well with our numerical results. 

It is generally accepted that the distribution of the de- 
cay rates 7 (at least for small 7) in the localised regime is 
log-normal, i.e., In 7 is distributed according to a Gaus- 
sian distribution. Recent interest has rather been in the, 
large-7 tail which was shown to follow a power-law.c3c3 
In a log- normal distribution, most of the weight lies in the 
right tail, so those papers give a sufficient description for 
most of the eigenmodes. In the context of applications 
to random lasers we are, however, interested in the small 
decay rate tail, hence in the log-normal distribution. 

To our knowledge, there is only a single paper by 
Titov and Fyodorov that gives explicit expressions for 
the parameters of that log-normal distribution.EJ How- 
ever, their analytical results are for a somewhat different 
system so it is difficult to tell whether they agree or dis- 
agree with our findings. We will return to this aspect 
at the end of this section. First, we want to present the 
results of our numerical simulations. 



Using the log-normal ansatz, the distribution of the 
decay rates 7 is 



P{"i) = foexp 



(Iog7-log7o) 



(13) 



The numerically computed histograms indeed follow this 
form, see Fig. ||, except for the large-7 tail — as already 
mentioned above but this deviation is only seen in a log- 
log plot. 

Fig. @ shows in the left the numerically computed 70 
as a function of L for N — lb. Also displayed is the lo- 
calisation length £ computed numerically from the trans- 
mission, being in good agreement with the analytical pre- 
diction (p^). The quality of the data is good enough to 
say with confidence that 70 decays exponentially with 
a length scale that is somewhat larger than £. Fig. (||) 
shows in the right the value of 70 also for two other values 
of N, and all three cases are well-described by introduc- 
ing a numerical fitting factor a, 

70 oc exp ( with a = 1.12 . (14) 

It is known that working at the centre of the conduction 
band when in the localised regime can introduce certain 
artefacts, especially in analytical approaches. Among 
other, the localisation length at the band centre can dif- 
fer by .-approximately 10% from the value outside the 
centrerj We have defined £ based on the transmission 
through the sample (at an energy corresponding to the 
band centre), and in transmission resonances at all ener- 
gies can contribute. A numerical prefactor a that differs 
by about 10 % from 1 thus does not come as a complete 
surprise. 

We still need to compute the proportionality factor ap- 
pearing in Eq. (|l4|). For this purpose we need to plot the 
ratio of the numerically computed 70 and the right-hand 
side of Eq. (|lj) for different values of N. We did this 
for L = 71.55 I. Since this is a very expensive operation, 
we have computed a large number of samples only for 
N = 10, 15, 20 so that their statistics is better than for 
the other values of N. An estimate of the error for these 
"better" data points has been included in the figure. This 
allows us to conclude that 



7o 



N 2 



exp 



— - ] with a 



1.12 



(15) 



It should be noted that this equation contains two nu- 
merical coefficients, and there is no obvious reason why 
they should be identical. Still, we find that they both are 
approximately a — 1.12. 

Re-introducing "physical units" into Eq. ( |l5| ) is a bit 
more difficult than it was for Eq. (§) where it was ob- 
vious that one simply has to multiply by the velocity of 
propagation c. Here one has to multiply by c/A where A 
is the perpendicular grid spacing. Due to the assumption 
of one propagating mode per (lateral) grid point made in 
Sec. [ilj A is not arbitrary but has a well-defined physical 
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FIG. 9: Left: Transmission through the sample and position 70 of the maximum of the log-normal distribution as a function 
of the length of the sample (for N = 15). There is a small but finite difference between the length scales for both quantities. 
Centre: Position 70 of the maximum of the log-normal distribution as a function of length for samples of different width N. 
Right: Prefactor in front of the exponential for L = 71.55 I for different N [cf. Eqs. @ and ^]. The dashed line marks the 
curve 1.12 /N 2 . 



meaning. For the electronic case, A = 7r/fcp with fep the 
wave vector at the Fermi level, and for the photonic case 
A = 2X/tt with A the wave length of the light (hence 
c/A = 1/(4*/)). 

Determining the width a of the distribution is more 
difficult since we can only use the left wing of the distri- 
bution — the right wing eventually turns into a power- 
law tail and thus no longer follows a log-normal distri- 
bution. Once again, we have accumulated more data for 
N = 10, 15, 20 so that some indication of the error is 
possible for those three data points. 

From our data, we propose the formula 



2/3 



(16) 



where a = 1.12 has the same value as in Eq. (fLq). As 
Fig. [l0| shows, there clearly is no disagreement between 
the numerical data and Eq. jl^). However, please re- 
member that the | should be thought of as a fitting fac- 
tor that might not be exactly 2/3 but perhaps rather 0.67 
or some other numerical factor. 

Since the distribution is log-normal only for not too 
large 7 (remember the power-law tail for large 7) the 
normalisation is nontrivial [^(7) is not normalised to 1 
any longer!] and cannot be computed from 70 and a. 
The constant a in Eq. ( |l3| ) is directly equal to the height 
of the peak of the numerically computed -P(/y). Since the 
total area underneath the numerically computed P(y) 
(and hence its normalisation) is dominated by the large-7 
tail, a has a relatively large error. Taking all the available 
data, the most likely value is 



(17) 



This value has been determined from a large number of 
simulations that for space reasons cannot be presented 
here. Unfortunately the quality of the data is not good 
enough to decide whether an additional prefactor a = 
1.12 should appear. 



3.5 
3.0 
2.5 
2.0 
1.5 
1.0 



1 1 1 


,, 1 1 1 ,,,,,, , 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i^j^i 
















7V = 20, 










: + 






:-F"' 




N= 10 '■ 




. 1 , , , . 1 . . . . 1 


, ... 1 .... 1 .... 1 .... 1 ... ,' 



35 40 45 50 55 60 65 70 75 

LI I 



M 0.66 




N 



FIG. 10: Top: Comparison of the numerically computed a 
with the prediction from Eq. (^). Bottom: Comparison of 
the numerically computed prefactor [from dividing the nu- 
merical a by (L/<) 2/3 , cf. Eq. @] with the prediction 2/3 
for samples of L — 71.55 I. For some values of TV a higher 
number of samples has been computed, so that the data qual- 
ity is high enough to also compute an error estimate. This 
has been indicated in the figure. 
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At the present it is not possible to tell whether our 
results agpe with the ones put forward by Titov and 
Fyodorov.li-3 In particular, they arrive at 



70 oc exp 



(18) 



whereas our finding ( |l5| ) was 70 oc exp(— L/1.12£). There 
are two obvious differences between the model used by 
them and the model employed by us. First, for numeri- 
cal reasons we work at the centre of the conduction band 
while they work near (but sufficiently far away from) the 
band edges. This might explain the factor a — 1.12 that 
we have to introduce. Secondly and probably more im- 
portantly, they consider a system of length L' that is 
closed at one end whereas our systems have length L and 
are open at both ends. It is obvious that a half-closed 
system of length L' corresponds to an open system of 
length L > L'. Eq. ( |l8| ) suggests that those two systems 
could be mapped into each other by setting L w 3L' but 
there is no further evidence to support this claim. 



VII. LASING THRESHOLD OF A RANDOM 
LASER 

A random laser is a laser where the necessary feed- 
back is not due to mirrors at the ends of the,4ase& but 
due to random scattering inside the medium.QCj'EJ We 
model the random laser as a disordered slab containing 
a dye that is able to amplify the radiation in a certain 
frequency interval with rate l/r a . The lasing threshold 
is the amplification rate at which the intensity of the 
emitted radiation diverges in a linear model. If satura- 
tion effects are included, the emitted intensity increases 
abruptly but finitely at crossing the lasing threshold. 

The lasing threshold is given by the value of the small- 
est decay, rate of all eigenmodes in the amplification 
window.Ea (Remember that 7 actually is twice the decay 
rate. On the other hand, also l/r a enters the relevant 
formulae only with a prefactor 1/2. 7 thus indeed gives 
the necessary amplification rate 1/V a .) This is easily un- 
derstood since this simply means that in the mode with 
the smallest decay rate the photons are created faster 
by amplification than they can leave the sample (=dc- 
cay). It, however, alaoJollows from a complete quantum 
mechanical analysis .@S 

The distribution of the decay rate has been computed 
in this paper. A certain number K of modes will be in the 
frequency window where amplification is possible. The 
lasing threshold is given by the smallest 7 of these K 
modes. In a simple picture that is valid once K 3> 1 we 
can assume that the K different 7's are distributed inde- 
pendently according to P(j)£3 The distribution P(7) of 
the smallest mode and hence of the lasing threshold then 
becomes 



P( 7 ) = KP( 7 ) 



1 - 



-1 K-l 



p(tW 



For K I, the distribution P(7) of the lasing thresh- 
old is not longer identical to the distribution P(7) of the 
decay rate of each individual mode. In particular, not 
only the precise form of these two distribution will be dif- 
ferent, but also the "typical" value of the lasing thresh- 
old can be different from the "typical" decay rate 70. 
Interestingly, for chaotic cavities in the diffusive regime 
it was found ithat the latter two quantities differ only 
insignificantlyLj'La which might seem counter-intuitive. 
A slab geometry is more "complicated" in that the scal- 
ing K oc N "tries" to lower the lasing threshold with 
increasing N. 

For K ^S> 1 the distribution P(7) is sharply peaked 
around its maximum. The position 7 m of the maximum 
is given by the solution of the equation dP(y m )/d'y m = 0, 
hence 







dP(jn 
d-Jm 



1 - 



-(K-l)[P( lia )f 



(20) 

While Eq. (|l^) is difficult to compute numerically due to 
the large exponent K — 1 3> 1, in Eq. (|2^) this exponent 
no longer appears. 

Eq. ( pfj| ) depends on P(7) which in turn depends on 
the dimensions L and N of the system. In assuming that 
the number of propagating modes is equal to the width 
N of the sample we already have made the assumption 
that the width (and hence also the length) is measured 
in units of A/2. (The "2" accounts for polarisation.) The 
total number of modes in the sample thus is LN. We as- 
sume that a fraction / of them is inside the amplification 
window of the dye, hence K = fNL. For simplicity we 
neglect complications as the shape of the mode profile. 
(It is easily incorporated into the numerics and we re- 
frain from doing this just to prevent having to introduce 
even more parameters.) / depends only on the chemical 
properties of the dye and not on the dimensions of the 
sample. 

In the following we will show how to compute the most 
likely lasing threshold for samples in both the diffusive 
and in the localised regime. 



A. Lasing threshold in the diffusive regime 

The change of the lasing threshold with increasing sys- 
tem size is influenced by a subtle interplay between L and 
N in determining the distribution P(7) and in determin- 
ing the number K = fNL of total modes. 

If K 3> 1 the lasing mode has a decay rate in the low-7 
tail of P(7) (i.e. 7 < 70 or y < 1). The weight of this 
tail is 



M 



M-l 

P(y)dy 
„ yy> y (M-l)! 



— M 



(21) 



(19) 



and goes to zero as M becomes larger. For M — > 00 the 
tail disappears completely, as is already obvious from the 
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asymptotic form of the distribution, 



P. 



M— >oo 



(y) = 



(y < 1) 

i/y 2 (i/>i) 



(22) 



With increasing A*" and hence increasing M, the prob- 
ability that a given mode has a small y thus decreases 
rapidly. On the other hand, we are interested in the 
smallest decay rate out of K modes, and K increases lin- 
early with N. This are two counter-acting effects, and it 
is not obvious which of these two is stronger. 

The effect of an increase of the system size L, on the 
contrary, is obvious. First, the average decay rate 70 
decreases according to Eq. (||). Secondly, M decreases 
from Eq. (|Io|), leading to even smaller values for 7 of the 
lasing mode. 

There have been some analytical atterants.to compute 
the lasing threshold for a chaotic cavitjOEj. For large 
M, the small- y tail of Eq. (0) was approximated by 



P(y) 



1 

2M 



l + crf(VM7%-l]) 



(23) 



This allows to arrive at scaling laws of the lasing thresh- 
old for variable M at fixed K. Unfortunately, the differ- 
ence between two counter-acting effects of an increase in 
N are so small that Eq. (^) is a bit too crude for our 
needs. 

We thus have to revert to a numerical procedure. 
Eq. (Q) can be rewritten using the incomplete Gamma 
function 



T(a, x) 



-*dt 



(24) 



r(a, 0) reduces to the well-known Gamma function T(a). 
For numerical reasons it is advisable to introduce the reg- 
ularised Gamma function Q(a,x) = T(a,x)/T(a). Fast 
numerical algorithms to compute Q(a,x) exist. [Please 
note that in the literature the definitions of the regu- 
larised Gamma function sometimes disagree in that our 
Q(a, x) is denoted as f — Q(a,x).] Now we can express 
Eq. (fjj) and its derivative and integral as 



p (y) = -3 [1 - Q( M + 1. M v)] > ( 25a ) 

dP(y) (My) M My 2, n(M 

(25b) 

/ P(y')dy' = - [Q(M + 1, My) - l] + 1 - Q(M, My) , 

Jo y 

(25c) 



B. Lasing threshold in the localised regime 
From Eq. (|I^) we directly arrive at 

(In 7 - ln 7 0) 2 



dP{l) , In 7 -In 70 

= — 2b - exp 



1 + erf 



(26a) 
2111(7/70) -<r 2 
2a 

(26b) 



A further simplification is not possible, and we did not 
manage to find suitable approximations. Also for the 
localised regime we thus are restricted to a numerical 
evaluation. 



C. Numerical results 

The lasing threshold is computed numerically from 
Eq. (p0|), using the formulae from Sees. VUA and V11B. 



Into the formulae presented there, we have to insert the 
correct dependence of the 70, M, er, etc. on L and A^ that 
was presented earlier in this paper. Despite this compli- 
cation the numerical calculation is straight forward as 
Eq. ( p0| ) possesses a single root only. Since this root has 
a change of sign, it is easily found numerically. 

Fig. [n] shows the results for both the diffusive and the 
localised regimes, for both f — 0.1 I and / = 0.001 I. 
(The mean-free path appears as a factor since the figure 
is in units L/l and not L.) The formulae found in this 
paper are valid deep in the diffusive regime respectively 
deep in the localised regime. Near the cross over, hence 
near the line L « Nl, this condition is not fulfilled. The 
numerical values near the diagonal line in Fig. [n] should 
thus be viewed with caution. 

As can be seen from the figure, in the diffusive regime 
with N ^ L/l the lasing threshold becomes almost in- 
dependent of the width N of the sample (for sufficiently 
large A^) , and the most likely value of the lasing threshold 
is about 



2d 
L 2 " 



(27) 



so that Eq. 



can be evaluated efficiently. 



hence the value given by Eq. (||). This means that even 
though K 1, P(y) for y < 1 is already so small that it 
dominates over the large value of K . Differences between 
this simple approximation and the precise numerical re- 
sult appear for finite N, with the size of this difference 
depending on /. However, for designing experiments it 
is obvious from the results presented here that the only 
feasible way to lower the lasing threshold of a random 
laser in the diffusive regime is increasing its length, not 
modifying its width. 

As Fig. [Tl] shows, also in the localised regime there is 
only a small dependence on /. This means that in a log- 
normal distribution the weight of the left tail is so small 
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FIG. 11: Most likely lasing threshold as a function of the length L/l and the width A^ of the sample. For L/l > N (the lower 
right part of the diagram) the sample is localised, for L/l < N (the upper left part) the sample is diffusive. The diagonal 
line marks the division, and the results close to that line should thus be viewed with caution. The left diagram depicts the 
results for / = 0.1, the central diagram for / = 0.001. The right diagram again depicts the results for the localised regime with 
/ = 0.1 but the horizontal axis has been rescaled to L/£ instead of L/l. The numbers x at the contour lines mean lQ x c/l in 
the diffusive regime, and lO^c/A in the localised regime. 



that unless K is exponentially large 7 m cannot become 
much smaller than the position 70 of the peak of the 
distribution. The difference to the diffusive regime is that 
the lasing threshold can be lowered efficiently not only by 
increasing the length but also decreasing the width N and 
hence driving the system farther into localisation. 

It is no surprise that samples in the localised regime 
generally have a lower lasing threshold than samples in 
the diffusive regime. We have shown that also the dif- 
fusive samples can have an "acceptably small" lasing 
threshold as it is trivial to make them very long (since 
there is no need to care much about their width). For 
both the diffusive and the localised regime, the typical 
decay rates of a single mode are comparable to the lasing 
threshold. 



VIII. CONCLUSIONS 

We have numerically computed the distributions of the 
residues (or decay rates) of a disordered slab. The slab 
has length L, mean free path I, width respectively cross- 
sectional area N (N is given as number of propagating 
channels) and velocity of propagation c. We were able 
to "guess" simple analytical formulae that are able to 
describe the numerical results well. 

For a sample in the diffusive regime (L < Nl) we found 



in Eqs 



■ a sample 
s. §|§ 



V 11 2lc V 2lc ' 
1 r r (! 

Hv) = 31- 



N Ny * 

1+L/Sl ' l+L/61 ' 



r(i 



JY 



l+L/61 ' 



(28a) 
(28b) 



where T(a,x) is the incomplete Gamma function. The 
agreement between the numerical results and the pro- 
posed formulae is good, and there is the possibility that 
Eq. ( pq ) might become exact in the limit L/l ^> N ^> 1. 
However, there is only numerical and no analytical evi- 
dence to back this claim. 

For a sample in the localised regime (L > Nl) with 
localisation length £ = (N + 1)1/2 we found in Sec. [v| 



L (Iog7-log7o) 2 



7o 



a 

N2 



at 


a 2 


) 








H) 













a = 1.12 , 



2/3 



(29) 



The quality of the simulations results in the localised 
regime is somewhat less than in the diffusive regime. For 
this reason, Eq. ( p9| ) should be understood as an approx- 
imate fit only, and it very probably differs from the exact 
relation, especially outside the band centre. 

These results can be applied to both electronic and 
photonic systems. For photonic systems we have shown 
that under realistic assumptions the lasing threshold of 
a random laser is close to 70 both in the diffusive and in 
the localised regime. Eqs. (£8|) and ( p9| ) thus not only 
give the distribution of the decay rate of each individual 
mode but also a good estimate of the lasing threshold, 
i. e. the smallest decay rate of a large number of modes. 
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